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Abstract 

The Vlasov-Poisson equations describe the evolution of a collisionless plasma, 
represented through a probability density function (PDF) that self-interacts via 
an electrostatic force. One of the main difficulties in numerically solving this 
system is the severe time-step restriction that arises from parts of the PDF asso- 
ciated with moderate-to-large velocities. The dominant approach in the plasma 
physics community for removing these time-step restrictions is the so-called 
particle-in-cell (PIC) method, which discretizes the distribution function into a 
set of macro-particles, while the electric field is represented on a mesh. Several 
alternatives to this approach exist, including fully Lagrangian, fully Eulerian, 
and so-called semi-Lagrangian methods. The focus of this work is the semi- 
Lagrangian approach, which begins with a grid-based Eulerian representation 
of both the PDF and the electric field, then evolves the PDF via Lagrangian 
dynamics, and finally projects this evolved field back onto the original Eulerian 
mesh. In particular, we develop in this work a method that discretizes the 1+1 
Vlasov-Poisson system via a high-order discontinuous Galerkin (DG) method in 
phase space, and an operator split, semi-Lagrangian method in time. Second- 
order accuracy in time is relatively easy to achieve via Strang operator splitting. 
With additional work, using higher-order splitting and a higher-order method 
of characteristics, we also demonstrate how to push this scheme to fourth-order 
accuracy in time. We show how to resolve all of the Lagrangian dynamics in 
such a way that mass is exactly conserved, positivity is maintained, and high- 
order accuracy is achieved. The Poisson equation is solved to high-order via 
the smallest stencil local discontinuous Galerkin (LDG) approach. We test the 
proposed scheme on several standard test cases. 
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1. Introduction 



The Vlasov equation in its various incarnations (e.g., Vlasov-Maxwell, Vlasov- 
Darwin, and Vlasov-Poisson) models the dynamics of collisionless plasma. Plasma 
is the state of matter where electrons have dissociated from their nuclei, creating 
a mixture of interacting charged particles. This mixture can evolve via a variety 
of effects, including electromagnetic interactions and through particle-particle 
collisions. In the collionlcss limit, the mean free-path is much larger than the 
characteristic length scale of the plasma; and therefore, particle-particle col- 
lisions are dropped from the mathematical model. Vlasov models are widely 
used in both astrophysical applications (e.g., [S1H1I37]), as well as in laboratory 
settings (e.g., HUHB M M HE] )■ 

The development of accurate and efficient numerical methods for the solution 
of the Vlasov equations are faced with a variety of numerical challenges, the most 
important of which we describe below. 

• High dimensionality. The Vlasov system is a nonlinear and nonlocal 
advection equation in six phase space dimensions (x 6 R 3 and v € R 3 ) 
and time - this of often referred to as 3 + 3 + 1 dimensions. Even though 
the Vlasov equation is in many ways mathematically simpler than fluid 
models, the fact that it lives in a space of twice the number of dimensions 
makes it computationally much more expensive to solve. 

• Conservation and positivity. In fluid models, conservation of mass, 
momentum, and energy are often relatively easy to guarantee in a numer- 
ical discretization, since each of these quantities is a dependent variable 
of the system. In Vlasov models it is generally more difficult to exactly 
maintain these quantities in the numerical discretization. Exact positivity 
of the probability density function is also not guaranteed by many stan- 
dard discretizations of the Vlasov system; and therefore, additional work 
in choosing the correct approximation spaces is often required. 

• Small time steps due to v e R 3 . In the non-relativistic case, the 
advection velocity of the density function in phase space depends linearly 
on the components of the velocity vector v € R 3 (see equation ^ in |2j). 
Since it is in general possible to have "particles" in the Vlasov system that 
travel arbitrarily fast, there will be a severe time-step restriction, relative 
to the dynamics of interest, that arises from parts of the PDF associated 
with moderate-to-large velocities. 

Several approaches have been introduced to try and solve some of these 
problems, including particle-in-cell methods, Lagrangian particle methods, and 
grid-based semi-Lagrangian methods. We briefly summarize each of these ap- 
proaches below. 

• Particle-in-cell methods. Particle-in-cell (PIC) methods are ubiquitous 
in both astrophysical (e.g., [40]) and laboratory plasma (e.g., [8]) applica- 
tion problems. The basic approach is outlined in the celebrated textbooks 
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of Birdsall and Langdon [7J and Hockney and Eastwood [57] , both of which 
appeared in the mid-to-late 1980s. Modern improvements to these meth- 
ods are still topics of current research (e.g., adaptive mesh refinement [40] . 
very high-order variants [2H 120], etc. . .). The basic idea is that the dis- 
tribution function is discretized into a set of macro-particles (Lagrangian 
representation), while the electromagnetic field is represented on a mesh 
(Eulerian representation). The main advantage of this approach is that 
positivity and mass conservation are essentially automatic, the small time 
step restriction is removed due to the fact that the particles are evolved 
in a Lagrangian framework, and the electromagnetic equations can be 
solved via standard mesh-based methods. The main disadvantages of this 
method are: (1) numerical errors are introduced due to the interpolations 
that must de done to exchange information between the particles and 
fields, and (2) error control is non-trivial since particles may either cluster 
or generate rarefied regions during the evolution of the plasma. 

• Lagrangian particle methods. One possible alternative to the PIC 
methodology is to go to a completely Lagrangian framework - this re- 
moves the need to interpolate between the particles and fields. Such ap- 
proaches are commonplace in several application areas such as many body 
dynamics in astrophysics [3J, vortex dynamics |31j . as well as in plasma 
physics [T2]. The key is that the potential (e.g., gravitational potential, 
streamfunction, or electric potential) is calculated by integrating the point 
charges represented by the Lagrangian particles against a Green's func- 
tion. Since the charges are point particles, evaluating this integral reduces 
to computing sums over the particles. Naive methods would need 0(N 2 ) 
floating point operations to evaluate all of these sums, where N is the 
number of particles, but fast summation methods such as treecode meth- 
ods [3_1 [3Tj and the fast multipole method [24] can be used to reduce this 
to 0(N log N). The main disadvantage of this approach is that it relies 
on having a Green's function, which for more complicated dynamics (i.e., 
full elect romagnet ism ) , may be difficult to obtain. 

• Semi-Lagrangian grid-based methods. Another alternative to PIC is 
to switch to a completely grid-based method. Such an approach allows for 
a variety of high-order spatial discretizations, and can be evolved forward 
in time via so-called semi- Lagrangian time-stepping. The basic idea is that 
the PDF sits initially on a grid; the PDF is then evolved forward in time 
using Lagrangian dynamics; and finally, the new PDF is projected back 
onto the original mesh. This gives many of the advantages of particle 
methods (i.e., no small time-step restrictions), but retains a nice grid 
structure for both the PDF and the fields, allowing extension to very high- 
order accuracy. There have been several contributions to this approach 
over the last few years. One of the first papers that developed a viable 
semi-Lagrangian method was put forward by Cheng and Knorr [9 . More 
recent activity on this approach includes the work of Parker and Hitchon 
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[55] , Sonnendriicker and his collaborators (see for example [301 HH GUI HH1 
EUSIEIISI), and Christlieb and Qiu (5J . 

The goal of the current work is to develop a high-order, grid-based, semi- 
Lagrangian method for solving the 1+1 Vlasov-Poisson equation. We focus on 
the 1+1 case, leaving the problem of high-dimensionality for future work. Our 
discretization is based on high-order discontinuous Galerkin representations and 
a high-order operator split semi-Lagrangian time-stepping method. We argue 
that this approach is a promising method that produces very accurate results 
at relatively low computational expense. 

This paper begins with a brief review of the Vlasov equations in f|2] Part of 
the focus of the present work is to develop a higher-order version of the classical 
Cheng and Knorr [5] operator splitting method, which we review in Sj3) The 
spatial discretization for the proposed method will be based on the discontinuous 
Galerkin method, which we briefly review in Sj4] The heart of this paper is [J5j 
which details all the aspects of the proposed method, and in particular, explains 
how to achieve high-order in space and time, mass conservation, and positivity 
of the distribution function. Finally in <|6] we apply the proposed scheme to 
a variety of standard test cases for the Vlasov-Poisson system, including the 
two-stream instability problem and Landau damping. 



2. Mathematical equations 

The Vlasov system describes the evolution of a probability density function 
(PDF) in phase space: 

f.(t,x,v) : K + x R d x U d -> R s , (1) 

where d is the spatial dimension and S represents the number of plasma species. 
This PDF denotes the probability of finding a particle of species s at time t, 
at location x, and with velocity v. Although the PDF is not itself a physical 
observable, its moments represent various physically observable quantities: 

p s (t, x) := / f s dv, (mass density of species s), (2) 

p s u B (t,x) :— / vf s dv, (momentum density of species s), (3) 

£ s (i,x):=-/ \\v\\ 2 f s dv, (energy density of species s) . (4) 



2 



Under the assumptions of a non-relativistic and collisionless plasma, the 
PDF for each species obeys the Vlasov equation, which is an advection equation 
in (x, v) phase space: 

|+vVJ s + ^(E + vxB). V v / S = 0. (5) 
at m s 
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The "particles" represented by this kinetic description do not interact through 
collisional processes; and instead, are only coupled indirectly through the elec- 
tromagnetic field. In general, the electromagnetic field satisfies Maxwell's equa- 
tions: 



(6) 
(7) 
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V • B = 0, V • E 



C 2 (7, 



where B is the magnetic field, E is the electric field, and the total charge density 
and total current densities are given by the following: 

q s T q s 



E 



-PsU s . 



(8) 



Note that the electromagnetic field variables, B(i,x) and E(t,x), as well as the 
mass and momentum densities, p s (t,x) and p s u s (t,x), only depend on time and 
the spatial coordinates x. 

In this work we will not consider the full Vlasov-Maxwell system for a many 
species plasma; and instead, we only consider the single-species Vlasov-Poisson 
equation. To arrive at the Vlasov-Poisson system, we start with Vlasov-Maxwell 
and assume that the charges are slow-moving in comparison to the speed of light; 
this allows us to replace the full electromagnetic equations with electrostatics. 
Furthermore, we consider only two-species: one dynamically evolving species, 
which we take without loss of generality to have positive charge and unit mass 
m = 1, and one stationary background species that has a charge of opposite sign 
to the dynamic species. Because the background charge is stationary, we will 
only need to solve a single-species Vlasov equation. These assumptions conspire 
to form the Vlasov-Poisson equations: 



f,t + v 
V 2 , 



/,« + • /, v = 0, (9) 
> = p(t,x)-p , (10) 
E, and — p is the stationary background 



where <j) is the electric potential: 
charge density. 

The Vlasov-Poisson system contains an infinite number of quantities that 
are conserved in time. These can be used as diagnostics in a numerical dis- 
cretization. We list four quantities that will be used in diagnosing our proposed 
scheme, all of which should remain constant in time: 



R d JR d 



\f\ dvdx, 



R d JR 



f 2 dvdx 



Total energy 
Entropy 



If f iwr/dvdx+l f i| E |i 2 dx, 

z JR d JR d z JR d 



/log(/) dvdx. 

R d JR d 



(11) 

(12) 
(13) 
(14) 
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Finally, we point out that in the current work we are concerned exclusively 
with the 1+1 dimensional version of the above equations with periodic boundary 
conditions in x. In this case, the Vlasov-Poisson system on £1 = (t,x,v) £ 
1R+ x [-L, L] x R is: 

f <t + vf, x + E(t,x)f, v = 0, (15) 
E <x = p(t,x) - p , (16) 

with periodic boundary conditions: 

f(t,-L,v) = f(t,L,v), E(t,-L) = E(t,L), and (j){t, —L) = <j)(t, L). 

In these expressions we used the shorthand notation: 

x := x , v := v 1 , and E := E . 

The total and background densities are 

/oo ^ />L 

f(t,x,v)dv and p := — J p(t,x)dx. (17) 

Note that po is m fad constant in time, due to conservation of mass on the 
periodic domain [— L,L]. 



3. Strang operator splitting 

If we momentarily freeze the electric field in time, the Vlasov equation (|9| 
can be viewed as an advection equation of the following form: 

/,* + a(v) • /, x + b(x) • / >v = 0. (18) 

Cheng and Knorr [5] realized that such an equation can be handled very effi- 
ciently if split into the following two sub-problems: 

Problem A: f, t + a(v) • /, x = 0, 
Problem B: f >t + b(x) ■ /, v = 0. 

The key benefit of this splitting is that each operator is now a constant coefficient 
advection equation (i.e., the transverse coordinate acts only as a parameter), 
each of which can be handled very simply with a variety of spatial discretization 
and semi-Lagrangian time-stepping. The down side of this approach, of course, 
is the introduction of splitting error. 

Cheng and Knorr [3] developed a second order accurate version of this scheme 
via Strang operator splitting 39J. Their scheme is summarized in Algorithm [lj 
It is worth pointing out that the electric field computed in Step 2, E n+ 5, is 
second order accurate in time, even though it is computed after advection in 
the x variables only. 
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Claim. Assuming that the current solution at time t = t n is known exactly, 
and that each step in Algorithm^is carried out exactly in space, velocity, and 
time, the density computed in Step 2 is second order accurate in time: 

p"+^p(V + ^,x)+0(At 2 ). 

This also implies that the electric field in Step 2 is second order accurate in 
time: 

E"+5 =e(V + ^,x) +0(At 2 ^ 



Proof. By assumption the PDF after the first step satisfies the following rela- 
tionship: 

/(x,v) ;= /(V,x-^v,v' 
We integrate this relationship in velocity to compute the density at time t n + ^ : 

:= J f (x, v) dv = J f (V, x - ^v, v) dv 
= J f (t n , x, v) dv - ^ V x • || v/ (t», x, v) dvj + 0(At 2 ) 
= p"-^V x -(p"u")+0(At 2 ). 
Finally, we use the fact that 

P n t = -V x (p n u») , 

in order to assert that 

= p» + + 0(At 2 ) = p (t + ^,x) + 0(At% 
which proves the claim. □ 



Algorithm 1 Cheng and Knorr 9 operator split algorithm. 

1. \At step on /, t + v/, x = 0. 

2. Solve V 2 </> = — p 0j and compute E n+ 2 = V0. 

3. At step on f tt + E"+3 ■ /. v = 0. 

4. I At step on /, t + v-/ x = 0. 
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The original method of Cheng and Knorr [9] employs a cubic spline spatial 
discretization. In the past few years, work on semi-Lagrangian solvers using 
Algorithm [l] with a variety of spatial discretizations, including modified cubic 
spline interpolants, has been carried out by Sonnendriicker and his collaborators 
(see for example [101 HH EH HH1 El HH [2H H] ) ■ Another recent contribution to this 
approach was the method of Christlieb and Qiu [3J] , who combined the operator 
split method with high-order WENO (weighted essentially non-oscillatory) finite 
differences for the spatial derivatives. 

The focus of the present work is to again consider operator splitting tech- 
niques, but this time with high-order discontinuous Galerkin (DG) method for 
the spatial discretization, and with fourth-order operator splitting techniques. 
We describe the basic DG framework in the next section <|4]and various details 
of both a second and a fourth-order in time operator split approach in |5] 



4. The discontinuous Galerkin (DG) method 

The modern form of the discontinuous Galerkin (DG) method was developed 
in a series of papers by Bernardo Cockburn, Chi- Wang Shu, and their collabo- 
rators [TBI US Ull ESI HI]- In this section we briefly review the DG method for 
a general two-dimensional conservation law on a Cartesian mesh. This section 
will also serve to introduce the notation that we will use throughout this paper. 

Consider a general 2D conservation law of the form: 

«,t + /(g,t,x), x + 5 (g,i,x) )V = 0, in xe!lcR 2 , (19) 

with appropriate initial and boundary conditions. In this equation q(t, x) 6 R m 
is the vector of conserved variables and f(q,t,x.),g(q,t,x) g R m are the flux 
functions in the x and y-directions, respectively. We assume that equation ( 19 1 
is hyperbolic, meaning that the family of m x m matrices defined by 

A <'- n >=-(i-i) T ' 20 > 

are diagonalizable with real eigenvalues for all x and q in the domain of interest 
and for all ||n|| = 1. 

We construct a Cartesian grid over fl = [a Xl b x ] x [a y , b y ], with uniform 
grid spacing Ax and Ay in each coordinate direction. The mesh elements are 
centered at the coordinates 

i ~ 2~) Ax and V] = ay + V ~ 2~) AlJ ' ^ 

On this grid we define the broken finite element space 

W h = {w h € ioo(fi) : w h \ T G P\ VT e %} , (22) 

where W h is shorthand notation for W Ax ' Ay . The above expression means 
that on each element T, w h will be a polynomial of degree at most q, and no 
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continuity is assumed across element edges. Each element can be mapped to 
the canonical element (£,17) G [— 1, 1] x [—1, 1] via the linear transformation: 

x = Xi + £—, y = Vj + r n—. (23) 

The normalized Legendre polynomials up to degree four on the canonical ele- 
ment can be written as 

^ = {l, V3Z, V3 V , 3ft,, ^(3£ 2 -l), ^(3ry 2 -l), 



»7(3r-i), -VW-i), V( 5 ? -30, V( 5f ? - 3? ?)> 



2 , v /, 2 „ v , n 2 ^ 2 



^(5^-30, ^(5^-3^ |(3^-l)(3r; 2 -l), 

105,4 _ 45 2 9 105 * _ 45, 9l 
8 ^ 4^ 8' 8' 4 ' 8 J ' 

These basis functions are orthonormal with respect to the following inner prod- 
uct: 

(v (m) , V (n) ) := j J J y m) &v)<P (n) &V)dt;dr l = 5 rnn . (24) 
We will look for approximate solutions of ( 19 ) that have the following form: 

M(M+l)/2 

t?) := Qv ) (t)<P (k Ht*v), (25) 

iJ fe=i 

where M is the desired order of accuracy in space (i.e., for fifth order: M = 5 
and M[M + l)/2 = 15). The Legendre coefficients of the initial conditions at 
t = are determined from the Z^-projection of q (x,y,0) onto the Legendre 
basis functions: 

Qlf(0):=(q h (0,^ V ),^ k Htv))- (26) 

In practice, these double integrals are evaluated using standard 2D Gaussian 
quadrature rules involving M 2 points. 

In order to det ermine the Legendre coeficients for t > 0, we multiply con- 
servation law ( 19 ) by the test function tpW and integrate over the grid cell 
7ij. After the appropriate integrations-by-part, we arrive at the following semi- 
discrete evolution equations: 

-0 W = C (t) (O t) ■— N w - 13 - (27) 
dt 3 13 W) ■ « Ax Ay ' 1 ' 
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where 



N 



AS. 



Arc 



-l J-i 

1 nf=l 
M tl„h 

f 

1 



e=-i 



.1}=! 



d£dr?, (28) 
(29) 
(30) 



rj=-l 



The integrals in ( 28 ) can be numerically approximated via standard 2D Gaussian 
quadrature rules involving (M — l) 2 points. The integrals in (29) and (30) can 
be approximated with standard ID Gauss quadrature rules involving M points. 



5. A high-order semi-Lagrangian DG method 

We describe in this section a semi-Lagrangian discontinuous Galcrkin method 
for solving the Vlasov-Poisson system. This method will have all of the following 
properties: 

1. Unconditionally stable; 

2. High-order accurate in space (5 th order); 

3. High-order accurate in time (4 th order); 

4. Mass conservative; and 

5. Positivity-preserving. 

All of these properties are explained in detail in this section. 

We begin by explaining the basic idea on the constant coefficient ID ad- 
vection equation in §5.1| The extension of this ID scheme to the 1+1 Vlasov 
equation via Strang operator splitting is described in section |5.2| A simple and 
efficient local discontinuous Galerkin solver for the Poisson equation is described 
in j]5.3| The generalization of the Strang-split scheme to higher-order splitting 



is shown in { 5.4 Basic properties including conservation of mass and positivity 



in the mean are proved in §5.5| Finally, in {5.6 a limiter that provides global 
pointwise positivity is described. 

5.1. A toy problem: the ID advection equation 

Consider the ID constant coefficient advection equation: 

f.t + vf, x = 0, (t, x) e R + x R, (31) 

with initial condition /(0, x). For simplicity of exposition, let's assume in the 
discussion below that that v > 0; the extension to the case v < is straight- 
forward. We consider solving this equation on a uniform mesh of elements, 7j, 
that each have width Aa;. We begin by projecting the initial condition onto the 
mesh: 

M 

^(0,0| =£^(0)^(0, (32) 
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(c) I 1 1 

Figure 1: Illustration of the shift + project method for solving the constant coefficient advec- 
tion equation in ID as described in j ]5.1| Panel (a) shows piecewise polynomial initial data; 
Panel (b) shows the initial data shifted by some amount (i.e., the exact evolution of the initial 
data); and finally, Panel (c) shows the solution after it has been re-projected back onto the 
original piecewise polynomial basis. 



where F h represents the finite dimensional approximation of f(t,x), M is the 
desired order of accuracy, and <£id(£) are the ID Legendre basis functions: 



¥>ffi={l. ^ ^(3£ 2 -l), | (35C 4 - 30C 2 + 3) 



A simple, high-order accurate, and unconditionally stable algorithm to update 
this solution can developed based on the following two steps: 

1. Exactly advect the initial condition over a time step At: 

f{t+At,x) = f(t,x-vAt) 

2. Project this solution back onto the mesh %. 

This process is illustrated in Figure[TJ These two steps can be compactly written 
for any starting time t n and final time t n+1 = t n + At as follows: 



^ +1 ) = ^£^U* n ) f 
fc=i J ~ i 

i M r 1 

1 k=l J-l+2u 



-l+2v 



^SU + 2-2i/)^(0de 



(33) 



where 



vAt 



and 



vAt 
Ax 



(34) 
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Here |_-J denotes the floor operation^] and < v < 1. By construction, update 
( 33 ) is unconditionally stable independent of the polynomial order of the spatial 



discretization. 



The integrals in equation (33) can be evaluated exactly. For example, in the 



case of piecewise constants and j = 0, (33 1 is nothing more than the first-order 
upwind scheme: 



F 



(l),n+l 



F; 



(1). 



),n 



(35) 



In the case of piecewise linear polynomials and j — 0, the scheme can be written 
as follows: 



R 



(l),n+l 



F (1) <" - v 



V3F, 



(2),n 



R 



(l),n 



v'sV (F, 



P (2),n 



R 



(2),n 



^(2),n+l _ F (2),n 



y/3u( RP' n -V3R 



n(2),n 



(2),n 



'3 j/ ( F 



F, 



(2),n 



2V3F>_;' 



2z/ 



(2),n 



F 



(2),n 



(36) 



(37) 



The above method is a close cousin to the Lax-Wendroff discontinuous Galerkin 
scheme (LxW-DG) of Qiu, Dumbser, and Shu [35]. In particular, in the case 
of piecewise linear polynomials the LxW-DG for the advection equation can be 
written as 



F 



F 



(l),n 



V3R 



(2),n 



R^ + VSR^I 



(2),n 



F 



(2),n+l 



3^2 (p(2),n _ p (2). 



F (2).n _ p (2),l 



F, 



(2),' 



(1), 



(2),Ti 



7 3F 



(2), 



F 



(1), 



7 3F 



(2),n 



(38) 



(39) 



We note that ( 38 )— ( 39 ) and (36)-(37) agree except in the v and v terms in 



the F^ 2 ) update. We argue below that these additional terms in the method 



given by (36)-(37) are crucial in ensuring stability up to CFL number one, and 



their absence in the LxW-DG method cause a non-optimal stability result. 



Claim. The numerical update given by (36 1 -(37 1 is stable for < v < 1 



Proof. We apply von Neumann stability analysis to (36)-(37) by assuming the 
following ansatz: 



2 this function takes a real input and rounds down to the largest integer that is smaller 
than or equal to the input. 
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where I = y— T. Plugging this into (36)-(37) yields 





n+1 


i?(2) 





i + f (C-i) 



V3i/(l-i/) (C-i) 



V3v(v - 1)(C - 1) 1 + 2^ 3 (1 - C) + 6^ 2 C - 3^(C + 1) 



i?(2) 



where C = e" 7 ^*. With some work, which is omitted here, one can show 
that the maximum modulus of the eigenvalues of the amplification matrix for 
< v < 1 is given by 



We note that 
which concludes the proof. 



g{v) — max( 1, 1 — Qv + 6z^ 
g{u) = l Vi/€[0,1], 



□ 



Claim. The numerical update given by (38 1 -(39 1 is stable for < v < J 



Proof. Using the same von Neumann ansatz as in the previous claim, this time 



applied to ( 38 )— ( 39 1, yields 





n+1 


i?(2) 





i + f (C-i) 

V3u(l - C) 1 



y/3u(l-v) (C-l) 
3i/(C + 1) + 3j/ 2 (C-1; 



jr(2) 



With some work, which is omitted here, one can show that the maximum mod- 
ulus of the eigenvalues of the amplification matrix for < v < 1 is given by 



ffLxW-DG(^) = max( 1, \-%v 

We note that 

5lx\v-dg(^) = 1 VV e 



but 3lxW-dg(^) > 1 Vi/ G ( ~,1 



which concludes the proof. Also note that (/lxW-dg^) an d g(y) are the same 
except for the additional term 6i/ 2 in g(y). Therefore, we conclude that this 
additional term is crucial in ensuring stability up to CFL number one for the 
method given by (36)-(37). □ 



In the above discussion we focused on the piecewise linear DG method. How- 
ever, more generally, the LxW-DG method behaves similarly to the standard 
Runge-Kutta DG methods (RK-DG) [16] in that the maximum allowable CFL 
number is inversely proportional to the polynomial order of the spatial dis- 



cretization. The modified LxW-DG scheme represented by ( 36 1—( 37 ) , and more 



generally by ( 33 1 , always has some modified and some additional terms in the 
update (e.g., the u 2 and v 3 terms in (37)) that produce a scheme with optimal 
stability. 
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(a) 



i-i-i.j.i 

Li-lJ.3 



T,,j 



(b) 




Figure 2: Illustration of the advection algorithm used in the proposed semi-Lagrangian 
method. Panel (a) illustrates a method with 3 Gaussian quadrature lines in the ^-direction. 
The solution along each line segment, e.g., Li i, in each element, e.g., Tij, is a ID polyno- 
mial. Each Gaussian quadrature line has a different velocity, and as such, each line will get 
shifted by a different amount; this is shown in Panel (b). The subscript labels A, B, and C in 
Panel (b) highlight the fact that after advection in the x-direction, each Gaussian quadrature 
line might come from a different elements. 



5.2. A semi-Lagrangian DG method for Vlasov-Poisson 

In this section, we describe in detail a semi-Lagrangian discontinuous Galerkin 
scheme for solving a quasi-lD problem of the form: 

f,t + a(v)f tX = 0, (i, x, v) e R + x Ex E. (40) 

This resulting scheme can then be inserted into Steps 1,3, and 4 of Algorithm 



for solving the 1+1 dimensional Vlasov-Poisson system given by (15) and (116 
We note that in Steps 1 and 4, a(v) = v; while in Step 3 the roles of x and v 
are reversed and a(x) = E n+ i(x). 

We construct a Cartesian grid over D, = [-L, L] x [— V mayi , V m a X ], with uni- 
form grid spacing Ace and Av in each coordinate direction. The mesh element 
Tij is centered at the coordinates 

Xi = -L + - Ax and v 3 = -V max + (j - An. 

each element can be mapped to the canonical element via the simple linear 
transformation: 

x = x t +£,—, v = Vj+r)—. (41) 

Next, we further subdivide each element by introducing for each horizontal 
row of elements, j, a set of M horizontal lines located at 

Vjk ■= Vj + m for h = 1, . . . , M, (42) 

where rjk are the roots of the M th degree Legendre polynomial. Along each of 
these lines, we pretend that we are solving a ID constant coefficient advection 
of the form: 

f,t + a(v jk )f tX = 0. (43) 
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1 . Forwards to locate interface 




2. Backwards to get solution at quadrature points 



Figure 3: Illustration of the forward and backward nature of the proposed semi-Lagrangian 
scheme. First, the cell edges are propagated forward from their initial time to their final time. 
Once these locations are known, Gauss-Legendre quadrature points are placed between the 
old cell edges and the new cell edges. In order to find solution values at these Gauss-Legendre 
points, we trace backwards along the characteristics to the initial time. 



This equation is then essentially solved by the ID method described in the 
previous section: §5.1[ This is depicted in Figure [2j where we have chosen 
M = 3 for illustration purpose^ Finally a fully 2D solution is reconstructed 
in each element, 7y, by summing the M ID solutions for k = 1, . . . , M with 
the appropriate Gaussian quadrature weights. This algorithm is summarized in 
Algorithm [5] 

In the semi-Lagrangian way of thinking there are generally two philosophies: 

1. Forward evolution: Each "particle"]^] is advected forward in time and then 
a solution at the the future time is reconstructed (e.g., Particle-in-cell 

and Crouseilles et al. [H]). 

2. Backward evolution: The solution at a specific location in the future is 
determined by tracing backwards along characteristics (e.g., Restelli et al. 

m)- 

We interpret the semi-Lagrangian DG method described in Algorithm [2] as a 
sort of mixed forward and backward method: 

1. Forward evolution phase: The cell edges are propagated forward from their 
initial time to their final time. This process determines the set of quadra- 
ture points necessary for the projection step. This forward evolution step 
is needed in order to attain mass conservation. 



3 The parameter M is chosen to coincide with the spatial order of accuracy of the scheme. 
In the numerical examples section we always take M = 5 to coincide with a fifth-order accurate 
discretization in x and v. 

4 A "particle" can be a quadrature point, element interface, etc. . . depending on the details 
of the particular method. 
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Algorithm 2 The proposed semi-Lagrangian algorithm for solving quasi-lD 
advection equations of the form f tt + a(v)f. x = 0. 
0. Start with the current solution: 



f h (t n ,x,v) 



M(M+l)/2 

E 



(44) 



where M is the desired order of accuracy in x and v. The canonical variables 



(£,77) are linearly related to the variables {x,v) via (41) 



1. In each element, construct M horizontal lines given by (42), where r\ k are 



the M Gauss-Legendre quadrature points. The solution in element 7y along 
each one of these lines is 



M(M+l)/2 



E 



(45) 



2. Advect the cell interfaces forward in time through ID advection along 
each Vjk- After forward advection along vj k , the i th cell will contain the old 
interface i — \— Ijk and this interface will be located a distance Aa; Vj k from 



the new interface i 



where 



ijk ■- 



a(vj k )At 
Ax 



and 



Vjk 



a{vj k )At 
Ax 



ijk- 



(46) 



This is the forwards phase illustrated in Figure [3j 

3. Next we trace characteristics backwards in time for each Vjk- The resulting 
process yields: 

M(M + l)/2 r -l+2u ]k 

m=1 1 (47) 

M(M+l)/2 



i V F (m) / 

A m=l J-l+2»j k 



where each of the above integrals are evaluated using ID Gauss-Legendre 
quadrature rules with M points. This step in the algorithm is the backwards 
phase illustrated in Figure [3] 

4. Finally, we update the solution by integrating in the vertical direction: 

M 

if )ineW = E^S' ( 48 ) 

k=l 

where ui k are the usual Gauss-Legendre quadrature weights for the M Gauss- 
Legendre quadrature points r\ k . 16 



2. Backward evolution phase: Once the old cell edge locations are known at 
the new time, Gauss-Legendre quadrature points are placed between the 
old cell edges and the new cell edges. In order to find solution values at 
these Gauss-Legendre points, we trace backwards along the characteristics 
to the initial time. 

This process is illustrated in Figure [3] 

5.3. Poisson solver 

We describe in this section how to efficiently solve the Poisson equation in 
Step 2 of the operator splitting approach shown in Algorithm [l] Consider first 
the ID Poisson equation oni£ [a, b] with mixed boundary conditions: 



p(x)-p , jX (a)=7, <j>(b)=0. 



(49) 



We apply to the Poisson equation the so-called local discontinuous Galcrkin 
method (LDG) (see [TJ [25] for two reviews of various approaches for solving 
Poisson equations via the DG method), and rewrite it as a system of two equa- 
tions: 



E 



P(x) - po, 
E(x). 



(50) 
(51) 



We expand 4>(x), E(x), and p(x) on each element as follows: 



E 



M 

{<t> h (x),E h (x),p h (x)- Po } =^r{#\ 

Ti k=l 

where M is the d esir ed ord er o f accuracy. 

We multiply (50) and (51) each by </4d(£) an< ^ integrate from £ 
i 1: 



Ax r iD \ 


i 

-i 


Aa; j 




i 

Ax 




i 

-i 


Ax 


/ w 

1 VlD.fi 



E h di 



-1 to 

(53) 
(54) 



Next, we apply the following one-sided rules in order to evaluate 4> h and E h at 
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the grid interfaces: 

M M 

4>\-l) := £ Vffi(-l) *< W = E(-!) fe+1 *i W - 

fc=l fc=i 

M A/ 

^(1) :=£ p<g(-l)*ffi =J2(-V k+1 V2k^l^ 

k=l k=l 

M M 



fc=l fe=l 

M M 

e\i) : = £ vSajif 5 - £ ^fc^r E 

k = l k=l 



(fc) 



Using these definitions, (53 1 and (54) can be rewritten as follows: 

M 



Y / V2k^lV2£-l (E| fe) + (-1) £ E^) - StkE[ k) = A, 

k=l 



k=l 



where S is an M x M matrix with entries given by 



Slk = J ^ID.C V>m d Z- 



M 



Note that the boundary conditions in ( 49 1 imply that 
E /l (a)= 7 = 

<f> h (b) = p =► ]>>i) fe+1 ^/2fc^T$, ( n fe > +1 =/3, 



fc=i 



,1/ 



fe=i 



where m x is the number of grid elements. 
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Putting everything together, (59) and (60) can be written in matrix form: 



1 

Ax 



1 

Ax 



A 

B A 
B A 



C D 
C D 



B A 



C 



where, for example, 



D 

C 



E 





Ei 






E 2 






E 3 






E ms _ 






' $! " 























p 3 



Ei 

E 2 



(64) 



E„ 



(65) 



E 



and A, B, C, and D are M x M matrices with entries given by: 



A ek = V2fc-lV2^-l - S ek , 
B tk = (-l)V2fc- l\/2^- 1, 
C ft = (-l) fe+£+ V2fc - In/21 -1-5, 



(-i) H V2fc - i\/2T^T. 



(66) 

(67) 
(68) 
(69) 
(70) 



The advantage of this formulation is that equations ( 64 ) and ( 65 ) are already 
in lower and upper triangular forms, respectively, and therefore can be easily 
be solved. The matrices A and C can be easily inverted once at the beginning 
of the calculation. 

5.3.1. Dirichlet boundary conditions 

The above method can easily be adapted to handle Dirichlet boundary con- 
ditions: 

<f>(a) = a, 4>{b)=P> (71) 

by noting that these boundary conditions are equivalent to the mixed BCs in 
(49) if we carefully choose the parameter 7 in (49). It can be shown that the 
correct choice for 7 is given by 



/3-a 
b — a 



i=l 



Ax 



poj ds 
Ax 2 



Vp (2) 

2^{b-a)jr[ 1 ' 



(72) 



where m x is the number of grid elements. 
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5.3.2. Periodic boundary conditions 

Periodic boundary conditions can also be readily handled: 



(a) = 0(6), E(a)=E(b), 



(73) 



by again noting that we need to carefully choose /3 and 7 in ( 49 1 . It can be 
shown that the correct choice for 7 is given by 



7: 



-i- / 8 ( p(s) - po) ds = Y x t vf ] + y Pf ) , (74) 



and (3 is arbitrary. Without loss of generality we simply take /? = 0. By solving 

1 the 

(75) 



(49) with 7 given by (74 1 and with j3 = 0, we obtain a solution (j>(x) with the 

(a) = 0(6) = 0. 



property that 



5.4- High-order split semi-Lagrangian method 

Now that the basic pieces are in place (i.e., a semi-Lagrangian solver for 
each split-piece ^5.2 and the Poisson-solver ^5.3), we are ready to introduce a 
fully fourth-order accurate method for the Vlasov-Poisson system. We describe 
in this section some important implementation details needed to achieve this. 
The resulting method is summarized in Algorithm [3j 



5.4-1. Fourth-order operator splitting 

Consider a time-dependent problem where the right-hand side is written as 
the sum of two differential operators A and B: 

q, t =A(q)+B(q). 

A fourth-order accurate operator splitting technique for such systems was devel- 
oped by Forest and Ruth [33] and by Yoshida [HI H2] • If we define the following 
two constants: 

71 = 1 /3 « 1.351207191959658, (76) 

2 i/3 

72 = - sa -1.702414383919315, (77) 
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then the fourth-order splitting approach of [231 S3 S3 can be written as a 
composition of the following seven stages: 

Stage 1: « 0.6756 At step on q, t =A(q), 

Stage 2: ^ At « 1.3512 At step on q, t = B(q), 

Stage 3: w ^ — ~ -0.1756 At step on q <t = A(q), 
Stage 4: 7 2 At « -1.7024 At step on g j4 = B(q), 

fry-, -l- ryr>)At 

Stage 5: ^ — — -0.1756 A* step on g jt = .A(g), 

Stage 6: 71 At « 1.3512 At step on q >t = B(q), 
Ti At 

Stage 7: w 0.6756 At step on q :t = A(q) . 

We note that this splitting approach requires some steps larger than At: Step 
2 and Step 6; as well as backward steps: Step 3, Step 4, and Step 5. 

5.4-2. Application to Vlasov-Poisson 

One difficulty with raising the temporal order of accuracy from two to four 
is the time-dependence of the electric field. In other words, the Cheng and 
Knorr [S] method does not completely reduce the Vlasov-Poisson to two constant 
coefficient problems, since the electric field remains time-dependent. In the case 
of Strang splitting, it turned out that one could easily generate a second order 
accurate representation of the electric field at the half time step, t™ + | At, as 
required in Step 3 of Algorithm [l] simply by carrying out Steps 1 and 2 of 
Algorithm [T] Additional attention must be paid in order to obtain temporally 
fourth-order accurate representations of the electric field. 

In order to avoid having to use the electric field at different points in time 
(i.e., a multi-step method), we construct the fourth-order Taylor polynomial 
centered at t = t n : 

E(t, x) := E" + (t - t") E» + - (t - t") 2 E» + i (t - t") 3 E» it . (78) 

The electric field value E" is computed from the Poisson equation: 

V • E" = VV l = P n - p . (79) 

The first time derivative of the electric field is proportional to the momentum: 

V-E 4 = p, t = -V-(pu) =*. E« = -(pu)", (80) 

and is thus readily computable. In order to compute the remaining time deriva- 
tives, we write down the evolution equations for the first few moments of the 
Vlasov-Poisson equation: 

p,t + V ■ (pu) = 0, (81) 
(pu) 4 +V-E = pE, (82) 
E, t +V-F = p(uE + Eu), (83) 
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where 



p := I fdv, pu := / v/dv, E := / vv/dv, and F := / vvv/dv 



Using equation ( 80 1 and the above moment evolution equations, we can compute 



the second and third time derivatives of the electric field entirely in terms of 
spatial derivatives: 



E,tt = -0ni) jt = V-E-pE, 
E )t « = V-E, t -p, t E-pE, t 

= V • (puE + pEu) - V • V • F • 



E V • (pu) 



2 

p U. 



(84) 
(85) 



It is clear from these expressions that in order to compute E itt and E^tt, we 
need to be able to compute first and second derivatives in space. One approach 
for doing this in the discontinuous Galerkin framework is to multiply by a test 
function and then integrate-by-parts. However, this approach will in general 
lead to a loss of accuracy. Instead, the approach taken in this work is to apply 
central finite differences that work directly on the Legendre coefficients of the 
function that needs to be differentiated. 

Consider the ^-projection of the function f(x), where / : ]R — > H, onto the 
space of piecewise polynomials of degree four on a uniform mesh of elements, 
71, that each have width Ax: 



/* =E*i ( %iS(0- 



(86) 



Therefore, f h represents the finite dimensional approximation of f(x). We can 
approximate the first and second derivatives of f(x) to O (Ax 5 ) accuracy by 
computing appropriate central finite differences of the Legendre cofficients . 
If we let D x f h and D xx f h represent the finite dimensional approximations of 
f'(x) and f"(x), respectively, then the central finite difference formulas on the 
Legendre coefficients are 



D X F, 
D X F 



D X F^ 



D X F„ 



D F 

^xx^i 

n c(5) 



2Ax 



- 2V5 Aiif 33 + 78 Axif^' 
Al if) - f V3V7A 1 F^ ) 
AiF^ - 14^5 AiF^ 
A lF ^ 



AiF 



1 

Ax* 



(5) 
(3) 



(87) 



A 2 F^> -V5A 2 F} 6) + 11A2F, 
A 2 F! 2) -IV3V7A 2 Fl 4) 
A 2 F (3) - 7V5A 2 F} 5) 



.(5)1 



A 2 F, 
A 2 F„ 



(4) 

>, 

(5) 
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Mesh 


f'(x) error 


log 2 (Ratio) 


f"(x) error 


log 2 (Ratio) 


25 


1.747 x 10~ 4 




8.292 x 10~ 5 




50 


5.543 x 10~ 6 


4.98 


2.672 x 10~ 6 


4.96 


100 


1.738 x 10~ 7 


5.00 


8.413 x 10~ 8 


4.99 


200 


5.437 x 10~ 9 


5.00 


2.634 x 10~ 9 


5.00 


400 


1.699 x 10" 10 


5.00 


8.364 x 10~ n 


4.98 



Table 1: Relative L2-norm errors for computing f'(x) and f"(x) using the Legendre coefficient 
finite difference formulas j87} and ( |88| , respectively. The example shown here is for f(x) = 
e sm(2nx) Qn a un ;f orm mesn D n < x < 1. Periodic boundary conditions are imposed to 

compute the derivative in the first and last elements: Fq = F^ and F^ +1 = F± for each 
fe = l,2,3,4,5. 



where 

A lF ^:=F^-F^, (89) 
A 2 F ( > k) := F$\ - 2Fl k) + f}% (90) 

To the best of our knowledge, this is the first time such formulas have been 
written down in the context of discontinuous Galerkin methods. In Table [I] we 
verify the order of accuracy by computing the first and second derivatives of 
f(x) — e sln ( 2 ^ x ) _ The errors in this table are computed using the relative 



errors defined by equation (B.4) with M = 5 and varying Ax. See Appendix B 
for more details. 

Once we have constructed an approximation to the time-dependent electric 
field, we are faced with an advection equation with time-dependent coefficients: 

/, t + E(i,x)./, v =0. (91) 

This equation can be readily solved to high-order via the method of character- 
istics. The key step in this approach is the evolution of the coordinates v as 
function of time: 

— = E(t, x) =>> v(t n + At) = v(t n ) + / E(t, x) (ft, (92) 

at J t n 

At 2 At 3 At 4 

v(f + At) - v(t n ) + At E" + — E™ + — E»„ + — E% t . (93) 



In other words, the semi-Lagrangian DG method as outlined in §5.2| remains 
largely unaltered by the fact that the electric is time dependent. The only 
difference is that interfaces and quadrature points are transported by equation 



( 93 ) instead of the simpler version of this equation when E is constant in time. 

We note than an important advantage of this approach is that, just as with 
Strang splitting, it requires only a single Poisson solver per time step. Finally, we 
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summarize the complete fourth-order splitting method for the Vlasov-Poisson 
system in Algorithm [3j 



Algorithm 3 Fourth-order operator split algorithm. 
1. Solve V 2 cj) = p n — Po and compute E" = V<; 



2. Compute 
E" = 



- 0u) n , 

V • E n - (pE) n , 



E n m = V • (puE + pEu)" — V • V • F" + E" V ■ (pu) n 



(p 2 »T 



where the spatial derivatives are computed via ( 87 1 and ( 88 ) 
Then construct 



E(i, x) := E™ + (t - t n ) E™ + i(t - t n f E« + i(i - i") 3 E« ti . 



3. ^Ai step on /, t + v • /, x = 0. 



4. 7l At step on / t + E (t, x) ■ /. v = 0, i G t" 



0, 7l At 



(71+72) 



At step on /t + v • /. x = 0. 



6. 72 At step on /, 4 + E (t, x) ■ /, v = 0, t G t n + 71 At, ( 7l + 72 )At 



(71+72) 



At step on /,t + v • /. x = 0. 



8. 7 i At step on /, t +E (t, x)-/, v = 0, t G t n + ( 7 i+72)At, (2 7l + 72 )At 



9. ^At step on /, t + v ■ f x = 0. 



5. 5. Mass conservation and positivity in the mean 

The description of the method so far has yielded a fourth-order accurate in 
time and fifth-order in space semi-Lagrangian method for the Vlasov-Poisson 
equations. It still remains to show that the method is mass conservative and 
that it is positivity-preserving. We prove mass conservation in the first theorem 
below. This is followed by a proof that each step in the operator split semi- 
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Lagrangian produces a solution that it is positive in the mearj^J 

Theorem 1 (Conservation). Each step of the semi-Lagrangian method described 
above is mass conservative. 

Proof. It suffices to show that the semi-Lagrangian scheme is conservative on 
the quasi-lD problem given by ( |40[ ) with periodic boundary conditions in x. 
Using the notation of Algorithm |2| the update for the mean value in cell Tij can 
be written in the form: 

M 

p (l),new _ V\ > <?W 



k=l 

M M(M+l)/2 r -l+2u jk 

= 2~E E ^ F i-l Ijkj ^ (m) (e + 2-2^ I?7fc )^ 

k=l m=l 

M M(M+l)/2 x 

+ ^E E "* F i-i kj <pM(Z-2»ik,Tlk)dZ. 
1 k=l m=l J-l+2m 

The total integral of f (t,x,v) over the entire computational domain can then 
be written as 

E^ ),now i £ UkF ™_ Ijh . r 1+2 "%M K+2 _ 2 ^ %)(lc 



l 

+ 2 



We make the following change of variables in the integrals above: 

s = £ + 2 — 2i/jj- and s = £ — 2i/jk, 
respectively, which yields: 



E p (l),now _ 1 , , p( m ) f 1 

i,j i,j,k,m 



i,j,k,m 



Since we are summing over all i, we shift the first index of F without changing 
the total sum; this step allows us to combine the two integrals into one: 



E^ ncw =2 E f 



3 We show in the next subsection how to turn positivity in the mean into global positivity. 
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Since (s , rj) is polynomial of degree at most M — 1 in 77, the Gaussian 
quadrature represented by the sum over k is exact: 

r-i pi 



^3 



ee^ } 



1 -j 



^ (m) (s,77) dsdr] 



= EE^ 

i.j m 

where 5i m is the Kronecker delta. 



□ 



Theorem 2 (Positivity in the mean) . Let M denote the spatial order of accuracy 
and let 

M 



K := 



(94) 



where [•] denotes the ceiling operatioi^ Let f h (t n ,x,v) be a function defined 
on the broken finite element space (122]) with q = M — 1, and let 



ti(t n ,Z,v) :=f h (t n ,x,v) 



(95) 



where (£,77) 6 [— 1,1] x [—1,1] are the variables on the canonical element. As- 
sume that fjj(t n , £, 77) is non-negative at all of the following 2MK points: 



&»0 = ($*,»?*) 

Vfc = 1,...,M and VI = 1, 



where $ jk := v jk (1 - s e ) + s e , 
Vjk (1 + st) - 1, 



(96) 
(97) 



, K . Ln the above expression, i>j k is given by 
(46), S£ denotes the I th quadrature point in the standard ID Gauss-Legendre 
rule with K points, and rjk denotes the k th quadrature point in the standard ID 
Gauss-Legendre rule with M points. 

If one time-step in one coordinate direction is taken using the semi-Lagrangian 
scheme as described in Algorithm^with f h (t n , x,v) as the initial condition, then 
the approximate solution at the end of this time-step will have a non-negative 
average in every element (independent of the time step At): 



F. 



(i). 



>o, v75,-en\ 



Proof. Using the notation of Algorithm [2j the update for the mean-value in 
element %j can be written as 



F. 



(l)v 



X E 
2 ^ 

k=l 
1 M 



-l+2u jk \ A/(M+l)/2 



m— 1 



F 



(m) 
-1-Ijk j 



^ m \i + 2-2v 3k ,r lk ) } dt 



k=l 



M(M+l)/2 

^ rf m l. <<P {m) (Z-2v jk ,r lk ) } d£. 



TCI — 1 



"this function takes a real input and rounds up to the smallest integer that is larger than 
or equal to the input. 
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We note that the terms inside the brackets are simply shifted solution values, 
allowing us to express the above update as follows: 

i M ( r-l+2vjk r 1 1 



where 



P>UO ■■= fli- Ijkj (t n ,Z + 2 - 2v jk , Vh ) for i G [-1, -1 + 2u jk ] , 
Pf jk (0--=fi- Iikj (t n ,Z-2v jk ,r 1k ) for £e [-1 + 2^,1]. 

Since P^ k {C) and P?j k {£) are polynomials of degree at most M—l, we can exactly 
evaluate each of the above integrals via Gauss-Legendre quadrature rules using 
K points (where K is defined in (94)): 

p(l),new 



F 

• J 



M ( K K \ 

= \ E E + E ^* (^*) 



where the tu^'s are the standard quadrature weights for Gauss-Legendre quadra- 
ture with K points. 

To conclude our proof, we note that since all of the quadrature weights in 
the above expression for F^p' ncw are strictly positive, we obtain positivity in 
the mean, 

new > n 
ij - U ' 



if /y(t",£j v) is non-negative at all of the 2MK points defined in (|96|-([97|. □ 

5.6. Positivity-preserving limiter 

One of the key assumptions in the above proof of positivity in the mean is the 
fact that solution prior to a time-step must be positive at all of points defined in 
(96)-(97). We show in this subsection how to the limit the solutions, including 
the initial condition, so that wc achieve positivity at all of these points. The 
key piece of technology necessary for achieving this positivity is a modification 
of the limiter of Zhang and Shu |44j . This limiter is simple to implement and is 
completely local to each element. 

The solution on some element T can be written as 

M(M+l)/2 

f\i,v)--= E F W <P W &V), (98) 

where M is the desired order of accuracy in space. We assume that the element 
average is non-negative: i 7 ^ 1 ) > 0. We sample this solution on a set of test 
points: 

(6, m) e [-1,1] x [-1,1] for i=l,2,...,P, (99) 
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and define: 

m:= min rfi). (100) 

2— 1,...,P 

Note that m G (—00, F^]. 

The limited solution is defined as follows: 

ry) := F« + ■ (/"(£, *?) - , (101) 

where 

# = min 1, — - . (102) 



m 



Note that < 9 < 1 and that 




if < m < F«, , . 

" ' 103 

if m<0, 

This means that if the solution is already non-negative at each of the test points, 
then this limiter does not alter the solution. On the other hand, if the solution on 
the element is negative at any of the test points, then the high-order corrections 
are damped until the solution is again non-negative. We are guaranteed that as 
9 — > 0, the solution will eventually become non-negative on the entire element 
since F (1) > 0. 

In practice we implement the limiting strategy as follows: 

1. During each of the stages labelled 3-9 in Algorithm^ we apply the positiv- 



ity limiter with the test points given by (96 1-(|97[). In this case P = 2MK, 
where M is the desired order of accuracy in space and K is defined by 
(94). As proved in Theorem [2j this guarantees that in each stage the 



approximate solution remains positive in the mean. 

2. After all of the stages of Algorithm [3] have been carried out, we apply 
the positivity limiter one more time to the solution, this time with the 
test points taken as the P = M 2 Gauss-Legendre quadrature points on 
[—1,1] x [—1,1]. This final limiting provides some additional positivity 
enforcement and allows us to compute a variety of integrals of the form 



(A.l) with function values that are non- negative. This is particularly 



entropy (A. 6 1 



useful in computing the ii-norm (A. 3 1, the total energy (A. 5), and the 



6. Numerical examples 

In this section we apply the proposed scheme to a variety of numerical test 
cases. We begin in §6.1| by considering a linear advection equation with a velocity 
field that produces solid body rotation. This example is primarily used to 
show the benefits of switching from second to fourth-order operator splitting 
strategies. In \Q.2 we verify the order of accuracy of the method on a forced 
Vlasov-Poisson equation where we know the exact solution. In the subsequent 



28 



three subsection we consider three standard problems for the Vlasov-Poisson 



system: ^ 6.3 the two-stream instability, §6.4| weak Landau damping, and S6.5 
strong Landau damping. Unless otherwise stated, all simulations below are done 
with 5 th order in space and with the positivity-preserving limiters as described 
in 95.61 turned on. 



6.1. Linear advection 

We first consider a linear advection under a divergence-free velocity field: 

g, t + u-Vg = 0. (104) 

We take the computational domain to be [0, 1] x [0, 1] and the velocity field to 
be solid body rotation around (0.5,0.5): 

u = (u(y),v(x)) = (tt (2y - 1) ,tt (1 - 2x)) T . (105) 

The initial condition is taken to be a smooth, compactly supported bump that 
is centered at (xo,yo) = (0.4,0.5): 

fcos 6 (5fr) if r<0.3, 
q(p,x,y) = l V3 > - ' (106) 

I otherwise, 

where 

r - y/(x-x ) 2 + (y-y ) 2 . (107) 

This problem, just as the Vlasov-Poisson system, is solved via operator splitting 
on the two operators: 

Problem A: q.t + u(y) q, x = 0, (108) 
Problem B: q. t + v(x) q tV = 0. (109) 

We run the initial condition out to time t = 1, at which point it should return 
to its initial state. The errors are computed using the relative L2 errors defined 



by equation ( C.4 1 with M = 5 and varying Ax = Ay. See Appendix C for 
more details. Convergence studies with Strang and the fourth-order operator 
splitting results are shown in Table [2j 

6.2. A forced problem: verifying order of accuracy 

Next we consider an example of a forced Vlasov-Poisson equation where we 
have an exact solution. The forced Vlasov-Poisson system is 

f,t + vf, x + El v =ip(t,x,v), (110) 

E iW = / f{t,x,v)dv-V^, (HI) 
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Mesh 


SL2 Error 


log 2 (Ratio) 


SL4 Error 


log 2 (Ratio) 


10 2 


3.215 x 1CT 1 




5.679 x 10- 1 




20 2 


7.185 x 1CT 2 


2.16 


3.113 x 10~ 2 


4.19 


40 2 


1.578 x 10~ 2 


2.19 


1.691 x 10~ 3 


4.20 


80 2 


3.923 x 1CT 3 


2.01 


1.010 x 10~ 4 


4.07 


160 2 


9.890 x 10~ 4 


1.99 


6.220 x 10~ 6 


4.02 


320 2 


2.454 x 10~ 4 


2.01 


3.843 x 10~ 7 


4.02 


640 2 


6.136 x 10~ 5 


2.00 


2.390 x 10~ 8 


4.01 



Table 2: Convergence study for the linear advection equation. Shown are the relative errors 
computed via ( |C.4[ | at time t = 1. All calculations were done with 5 th order accuracy in 
space using the positivity-preserving limiters and a CFL number of 5.00, where CFL := 
Almax{max 9 \u(y)\/ Ax, maxj \v(x)\/Ay}. SL2 refers to the Strang split semi-Lagrangian 
scheme and SL4 to the fourth-order split semi-Lagrangian scheme. 



on (t, x, v) € [0, oo) x [— 7r,7r] x (—00,00) with periodic boundary conditions in 
x. We take the following source term: 

if)(t,x,v) = 1 sin(2:r-27rf)e~i (4t '- 1)2 ( (20F+1) Uv - 2^) 

2 I (n2) 

- ^/tt(Av - l)cos(2a; - 2nt) L 

The exact solution in this case is 

f(t,x,v) = {2-cos(2a;-27rf)}e-j( zk '- 1 ) 2 , (113) 

E(t,x) = sin (2a; - 2irt) . (114) 

The numerical scheme for this forced problem is the same as the one de- 
scribed in ^J5] with two minor modifications. First, the two operators in the 
operator split scheme are 

Problem A: f, t + v f, x = i/>(t, x, v), (115) 
Problem B: f tt + E(t, x) f tV = 0, (116) 

which means that Problem A is slightly modified from the unforced Vlasov- 
Poisson system. The modified A still has the same characteristics as the case 
with no source term; the only difference is that the solution is no longer constant 
along the characteristics. In order to advance / forward under the influence 
of operator A, we use the method of characteristics and obtain the following 
solution: 

f{t n+1 ,x,v) = f(t n ,x-vAt,v) + J ^(s,x + v(s-t n+1 ),v)ds. (117) 
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The time integral in the above expression can be easily exactly evaluated; we 
omit the details here. 

The second modification comes from the fact that with a non-zero source 
term, it is no longer true that E t = —pu. This means that the time interpolation 



described in { 5.4.2 must be slightly modified. Instead of using E t = —pu, (84) 



and (851, we make use of the following modified formulas: 



E tt (t,x) = -pu + C 1 , (118) 
E, tt {t,x)=E, x -pE + C 2 , (119) 
E, m (t, x) = (2puE) x - F ^ + E {pu)^ x + p 2 u + C 3 , (120) 



where 



Ci := — + — (4tt- l)cos(2a; - 2nt), (121) 
4 8 

C 2 := 3 ^ + in i6 16V ^ sin(-2x + 2nt) + ~ sin(4x - Ant), (122) 

„ 7T 7\fn + 167T — 64V7T^ . „ . 37T . , . ,,„„ N 

C 3 := -- + — — cos (2x - 2nt) - — cos 4x - Ant). 123 

4 32 16 

In the above expression we used the shorthand notation: 

x — x 1 , E:=E X , u:=u l , E := E 11 , and F:=F m . 

We run the initial condition out to time t = 1, at which point it should 
return to its initial state. Convergence studies on various grids on the domain 
(x, v) € [— 7r,7r] x [— 7r, 7r] with Strang and the fourth-order operator splitting 
results are shown in Table [3] We compute the errors in an identical manner to 
the linear test problem presented in the previous section using the relative Li 
errors defined by equation ( |C.4[ ) with M — 5. See Appendix C for more details. 



6.3. Two-stream instability 

The two-stream instability problem has become a standard benchmark to 
test numerical Vlasov solvers, and has been used as such by several authors 
(e.g., [ISlUlEllilllllH])- We use the following initial distribution function 

f(t = 0, x, v) = -j= (2 - cos (I)) e-£ , (124) 



8tt 

and solve on the domain (x, v) G [— 2ir, 2tt] x [— 2ir, 2tt]. Results for time t = 45 
are presented in Figure [4] for various mesh sizes. In Figure [5] we present vertical 
cross-sections of the solution taken at x = for the same mesh sizes. 

Figures |4] and [5] clearly demonstrate the ability of the discontinuous Galerkin 
methodology to approximate very rough data, something that is more difficult 
with methods that act over larger stencils. The results shown in these figures 
indicate far more structure than what is shown in other recent work, including 
results from the WENO method [34], [2] and an explicit DG method that uses a 
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Mesh 


SL2 Error 


log 2 (Ratio) 


SL4 Error 


log 2 (Ratio) 


10 2 


5.210 x 10- 1 


- 


9.493 x lO" 1 


- 


20 2 


1.433 x 10" 1 


1.86 


2.715 x 10" 1 


1.81 


40 2 


1.640 x 10~ 2 


3.13 


1.652 x 10~ 2 


4.04 


80 2 


3.438 x 10~ 3 


2.26 


7.079 x 10~ 4 


4.55 


160 2 


8.333 x 10~ 4 


2.04 


3.434 x 10~ 5 


4.37 


320 2 


2.068 x 10~ 4 


2.01 


1.962 x 10~ 6 


4.13 


640 2 


5.161 x 10~ 5 


2.00 


1.203 x 10~ 7 


4.03 


1280 2 


1.290 x 10" 5 


2.00 


7.509 x 10-° 


4.00 



Table 3: Convergence study for the forced Vlasov-Poisson equation. All calculations presented 
here are 5 th or der in space and were run with a CFL number of 5. Shown are the relative errors 
computed via §CA\ at time t = 1. SL2 refers to the Strang split semi-Lagrangian scheme and 
SL4 to the fourth-order split semi-Lagrangian scheme. Since the positivity-preserving limiters 
as described in §5. 6| don't guarantee positivity in the mean in the presence of a source term, 
we have turned them off for this convergence study only. 



piecewise constant representation of the distribution function, /, and a piecewise 
quadratic representation of the electric potential <f) |25| . 

In Figure [6] we demonstrate the effects of adding the posivitity-preserving 
limiter. We see that even without limiting, the base scheme already does a rea- 
sonable job of not producing large negative values in the distribution function. 
With the positivity-preserving limiters we are able to remove these small posi- 
tivity violations. In Figure [7] we plot four quantities that are exactly conserved 
by the continuous Vlasov-Poisson equation, but only approximately conserved 
in our numerical discretization: Li-norm of / (11), L2-norm of / (12 1, total 
energy (13), and total entropy (14). In particular, we use the numerical ap- 



proximations to ( 11 )— ( 14 ) as given by equations (A.3)-(A.6) in Appendix A 



We note that it is difficult to obtain accurate values for the total entropy (14), 
because there are many values where / becomes very small. 



6.4- Weak Landau damping 

Landau damping has been extensively studied both numerically US] 
and analytically (e.g., the work of Mouhot and Villani (32]). Just as the two- 
stream instability problem, Landau damping has become a favorite standard 
test case. It is particularly useful since the linear decay rates of the L 2 -norm of 
the electric field are well-known. 

We use the following initial distribution function 



f(t = 0,x,v) = 




(125) 
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f(1,x,v) at I = 45 



f(1,x,v) at t - 45 





f(1,x,v) at t = 45 




f(1,x,v) at t = 45 



0.35 
0.3 
0.25 
0.2 
0.15 
0.1 
0.05 





I 



0.35 
0.3 
0.25 
0.2 
0.15 
0.1 
0.05 




f(1,x,v) at I = 45 



f(1,x,v) at I = 45 





1 



0.35 
0.3 
0.25 
0.2 
0.15 
0.1 
0.05 




Figure 4: The two-stream instability problem. The panels in the left-hand column are results 
using the 2 nd order Strang splitting method. The panels in the right-hand column are results 
using the 4 th order splitting method. All simulations are 5 th order in space. The mesh 
sizes for the first, second and third rows are (m x ,ra v ) = (65,65), (m x ,rn v ) = (129,129), and 
(m x ,m v ) = (255, 255), respectively. All solutions use the positivity-preserving algorithm. The 
above figures were produced by plotting the numerical solution at each of the 5x5 Gaussian 
quadrature points in each mesh element. 



33 



f(t, x = O.OOO.v) at t = 45 



f(t, x = O.OOO.v) at t = 45 




f(t, x = O.OOO.v) at t = 45 



f(t, x= O.OOO.v) at t = 45 






-6 


-4 


-2 2 
f(t, x= O.OOO.v) at t = 45 


4 


6 



f(t, x= O.OOO.v) at t = 45 




Figure 5: The two-stream instability problem. The panels in the left-hand column are results 
using the 2 nd order Strang splitting method. The panels in the right-hand column are results 
using the 4 th order splitting method. All simulations are 5 th order accurate in space. The 
mesh sizes for the first, second and third rows are (m x ,m v ) = (65, 65), (m x , m v ) = (129, 129), 
and (m x ,m v ) = (255, 255), respectively. All solutions use the positivity-preserving algorithm. 
Each figure above represents the numerical solution at x = and use 5 Gaussian quadrature 
points for each cell in the D-coordinate. The above solutions were computed with an odd 
number of elements in each coordinate direction in order to easily obtain a slice of the solution 
at x = 0. 
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f(t, x = 0.O00,v) at t = 45 



f(t, x= Q.O00,v) at t = 45 




Figure 6: The two-stream instability problem. These panels demonstrate the effect of the 
positivity preserving-limiter. The result on the left is the 4 th order splitting algorithm with 
limitcrs turned on, and the result on the right hand side is the same algorithm with the 
limiters turned off. Both pictures represent a slice of the solution at x = and final time 
t = 45. Both results represent a mesh of size (m x ,m v ) = (129, 129) and are 5 th order accurate 
in space; an odd number of grid elements are used in order to easily obtain function values. 
We further note that mini /(^ii^i) = —2.020 X 10~ 2 for the solution without the limiter and 
mini f(%i,Vi) = 7.000 X 10 — 12 for the limited solution, where the minimum is taken over all 
5x5 Gaussian quadrature points over every mesh element. 



with a = 0.01 and k — 0.5 on the domain (x, v) E [—2ir, 2ir] x [— 2tt, 2w}. Because 
a is a small parameter, we expect to see results that closely agree with the linear 
theory, where the electric field decays exponentially. In Figure [8] we present this 
decay provided by 




log(||£(V)lk) :=logW / \E{t,x)Y dx 



versus time for two different mesh sizes. Our computed decay rate matches the 
linear decay rate, 7 = —0.1533. In Figure [9] we again plot the deviations of 
several quantities that are conserved by the continuous Vlasov-Poisson system 
from their initial values: H/Hij, ||/||l 2 i total energy, and entropy. We again use 



the numerical approximation to these norms given by (A.3|-(A.6) in Appendix 
[A"| Our results are comparable to what is reported for example by Qiu and 
Christlieb [34]. 

6.5. Strong Landau damping 



The initial condition is again given by (1251, this time with a = 0.5 and 
k = 0.5 on the domain (x, v) e [— 2ir,2n] x [— 2n, 2ir]. The time evolution of 
the distribution function is shown in the panels of Figure |10| These images are 
comparable to what is shown in Qiu and Christlieb |34| . but we are again able 
to capture more fine scale structure with the discontinuous Galerkin approach. 

A semi-log plot of the L2-norm of the electric field is provided in Figure 
|11| and decay rates are computed by sampling the solution at data points. 
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' 1 5 10 15 20 25 30 35 40 45 °' 1 5 10 15 20 25 30 35 40 45 



Figure 7: The two-stream instability problem. Shown in these panels are the deviations of 
the L\ norm of / (top-left), L2 norm of / (top-right), total energy (bottom-left), and total 
entropy (bottom-right) from their initial values. All simulations use a constant CFL number 
of 2.0, and are obtained from a mesh of size (m x ,m v ) = (129,129). The domain for this 
problem is (x,v) £ [—2-k,2tt] x [— 2ir, 2tt]. Each simulation is 5 th order accurate in space and 
is positivity preserving. 



3G 



^(Electric Field) ^(Electric Field) 




Figure 8: The weak Landau damping problem. Shown in the panels are semi-log plots of the 
Z/2-norm of the electric field. Simulations use a constant CFL number of 2.0 and are 5 th order 
accurate in space. All simulations use the positivity-preserving limiter. The figure on the left 
represents a mesh of size (m x ,m v ) = (64, f 28) and the result on the right was represents a 
mesh of size {m x ,m v ) = (f28,256), both on the domain (x,v) £ [— 27r,2-7r] X [—2tt,2tt]. Both 
simulations match the theoretical decay rate of 7 = —0.1533, and to demonstrate this we plot 
the line defined by y = 0.06e 7 '. One should note that the discrepancy in the two plots is due 
to the fact that twice as many time points in the second plot as the first one. 



We find that the initial linear decay rate is approximately 71 ps —0.292 which 
is close to the value of —0.243 computed by Zaki et al. [43], closer still to 
the value of —0.281 computed by Cheng and Knorr [5], but much larger than 
the value of —0.126 computed by Heath et. al 25J. In this same figure we also 
estimate the growth rate due to particle trapping and find it to be approximately 
72 = 0.0815; this number also differs from the value reported by Heath et al. [23] : 
72 = 0.0324. The initial linear decay was computed by taking the maximum of 
the first two peaks located at t « 2.45 and t 4.54. For the particle trapping 
growth regime, we sampled the maximum of the solution at the two peaks 
located at t w 2.33 and t rs 2.84. We postulate that the difference between 
our computed growth rates and those of Heath et al. J25] stems from the fact 
that we are using piecewise quartic polynomials to represent the distribution 
function, while they are using only piecewise constants. This issue should be 
further investigated. 

In Figure [12] we again plot the deviations of several quantities that are 
conserved by the continuous Vlasov-Poisson system from their initial values: 
ll/IUi, ||/||l 2 j total energy, and entropy. In particular, we use the numerical 



A.3 


HA.6 


) in 


Appendix 



|A| Our results are comparable to what is reported for example by Qiu and 
Christlieb [31]. 

Finally, we note that in both our weak and strong Landau damping simula- 
tions we made V max = 2tt instead of the more commonly used value of Vm ax = 5. 
The reason for this is that we noticed that between roughly v — 5 and v = 6 
the distribution function f(t,x,v) still had a non-negligable amplitude on the 
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Figure 9: The weak Landau damping problem. Shown in the panels are the deviations of 
various quantities that are conserved by the exact equations from their initial values. All 
simulations use a constant CFL number of 2.0, are 5 th order accurace in space, and are 
positivity preserving. The mesh size for the left hand column is (m x ,m v ) = (64,128) and 
the mesh size for the right hand column is (m x ,m v ) = (128,256). We note that the largest 
deviation for total energy for the 4 th order algorithm is on the order of 10~ 10 . 



order of about 10 6 ; the precise behavior of strong and weak Landau damp- 
Therefore, truncating at V m 



ing in this region is shown in Figure 13 Therefore, truncating at V ma , x = 5 
caused additional errors when tracking various conserved quantities; we found 
improvements in these errors when taking V max — 2tt. 



7. Conclusions and future work 

We have described in this work a semi-Lagrangian discontinuous Galerkin 
method for solving the 1+1 Vlasov-Poisson system. This method was shown to 
have all of the following properties: 

1. Unconditionally stable; 

2. High-order accurate in space (5 th order); 

3. High-order accurate in time (4 th order); 

4. Mass conservative; and 

5. Positivity-preserving. 

The proposed method is based on a series quasi- ID semi-Lagrangian advec- 
tion steps coupled with a fourth-order operator splitting scheme. The spatial 
discretization is handled via high-order discontinuous Galerkin representations. 
The Poisson equation is solved to high-order via a modified local DG method, 
where the boundary conditions are set so that the discrete Laplacian matrix 
is by construction LU factored. The scheme was applied to several standard 
Vlasov-Poisson test cases, which demonstrated the accuracy and robustness of 
the proposed scheme. 

The advantage of DG over other methods that are based on larger stencils is 
its ability to represent very rough data. We showed that this feature is important 
in the case of the two-stream instability and the Landau damping calculations 
presented in ^j6] With an explicit time-stepping method, the price that is paid 
for this spatial localization is a maximum CFL number that decreases with the 
spatial order of accuracy. In the context of Vlasov-Poisson we have tamed this 
problem by using the semi-Lagrangian framework. 

Future work will focus on extending the results described in this paper to 
higher-dimensional Vlasov-Poisson equations. Furthermore, modifications of the 
current approach to both the non-relati visit c and the relativistic Vlasov-Maxwell 
equations will be considered. 

Acknowledgements. The authors would like to thank the anonymous review- 
ers for their valuable feedback. This work was supported in part by NSF grants 
DMS-0711885 and DMS-1016202. 



Appendix A. Numerical evaluation of conserved quantities 



The conserved quantities defined in ( [11] ) (ii-norm), ( 12 ) (L2-norm), ( 13 ) (to- 
tal energy), and ( 14 ) (entropy) are used as diagnostics of the numerical methods 
proposed in this work. 
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Figure 10: The strong Landau damping problem. Shown in the panels are the distribution 
function at various points in time. This simulation was run with a constant CFL number 
of 2.0 on a mesh of size (m x ,m v ) = (128,256) using 5 th order accuracy in space and the 
positivity-preserving limiters. It is clear from these plots that the high-order discontinuous 
Galerkin method is able to capture much of the fine-scale structure for the solution. 
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L (Electric Field) ^=-0.2952, y z - 0.0844 L 2 (Electric Field) y = -0.2920, y 2 =0.0815 




Figure 11: The strong Landau damping problem. Shown in these panels are semi-log plots of 
the Li norm of the electric field with two different resolutions; the mesh size for the the figure 
on the left is (m x , m„) = (64, 128) and the figure on the right is (m x , m v ) = (128, 256). Both 
simulations use the positivity preserving limiter, are 5 th order accuracy in space and use a 
constant CFL number of 2.0. In each panel, 71 refers to the slope of the initial decay, and 72 
refers to the growth rate between times t = 20 and t = 40. 



In order to evaluate all of these conserved quantities in the numerical evo- 
lution, wc define the following functional: 



. . m x m„ M 2 

i h ( 9 (f h )) -^EEE^f/* 



i=l j=l k=l 



& Ax 



T] k Av 



(A.l) 



where m x is the number of elements in the x-direction, m v is the number of 
elements in the ^-direction, and ui). and are the M 2 Gauss-Legendre 

quadrature weights and points on [—1,1] X [—1,1], respectively. Expression 
(A.l) gives a numerical approximation to integrals of the form: 



I (9(f)) ■= 



L poo 



9(f(x,v)) dvdx. 



(A.2) 



L J —00 



Using (A.l) we define the following numerical approximations to the norms 

(A.3) 



dehned by ( 11 )— ( 14): 



11/1 



11/1 



L 2 • = 



Ih (\f h \) 

Ax Av 

~4~ 



m x m„ M(M+l)/2 

EE E 

i=i j=i i = \ 



M 



Total energy 
Entropy 



E 



(0 



i=l 1=1 



-I h (f h log(/' 1 )) . 



(A.4) 

(A.5) 
(A.6) 
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Figure 12: The strong Landau damping problem. Simulation results for the L\ norm (first 
row), Z/2 norm (second row), energy (third row), and entropy (bottom row) for strong Landau 
damping. All simulations use a constant CFL number of 2.0 and are 5 th order accurate in 
space. The mesh size for the left column is (m x ,m v ) = (64, 128) and the mesh size for the 
right column is (m x ,m v ) = (128,256). 
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Figure 13: A comparison of vertical slices for strong landau damping (left) and weak Landau 
damping (right) at time t = 60. The simulations for each panel use a grid resolution of 
(m x ,m v ) = (128,256), each are 5 th order accuracy in space, and each use the positivity- 
preserving limiter. The CFL number for each simulation is 2.0. We note that the solution 
is non-zero for \v\ > 5 in both cases, although there is much more activity in the case of 
strong Landau damping. The plots suggests that the commonly used maximum velocity of 
| p| = 5 should be increased in order to get better accuracy in conservation of the quantities 
JA.3HA.6b. 



Appendix B. Relative Z/2-norm error in ID 

Let f(x) be the exact solution of some problem of interest. Let f h (x) denote 
an approximation to f(x) using a discontinuous Galerkin method. On each 
element f h (x) and f(x) can be written as 



M 



Ti k=l 
Ti k = l 



(B.l) 



(B.2) 



respectively. The relative L2-norm of the difference on the domain x G [a, b] 
between the approximation, f h (x), and the exact solution, f(x), is given by 



\f(x)-f h (x)\\ 



J b a [f(x)-f h (x)] 2 dx 
J b a f(xrdx 



F (k) _ j.{k) 



2 \ 2 
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where N is the total number of grid elements and Ax = (b—a) /N . Therefore, we 
take as our relative L2-norm indicator the following easily computable quantity: 



E 2 (Ax,M) :-- 



EN 
»=1 l^k=\ 



?(*>) 



(fe)l 



EN 



fc=l 
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U 1 
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Appendix C. Relative Z/2-norm error in 2D 

Let f(x,y) be the exact solution of some problem of interest. Let f h (x,y) 
denote an approximation to f(x, y) using a discontinuous Galerkin method. On 
each element f h (x,y) and f(x,y) can be written as 



f h (x,y) 
f(x,y) 



M(M+l)/2 

= E *;?V fe) (^)> 

oo 

Tii fc=l 



(C.l) 
(C.2) 



respectively. The relative L 2 -norm of the difference on the domain (x, y) e 
[a x ,6 x ] x [a y ,b y ] between the approximation, f h (x,y), and the exact solution, 
f{x,y), is given by 



||/(s,y)-/ h (s,y)|| La _ f/a:/«S'[/( a; 'I')-/ h ( a: »»)] 2d » da: 



ll/O^lh 
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} + 0(Ax M ,A 2/ M ) 



where N x and iV^ are the the number of grid elements in each coordinate direc- 
tion, Ax = (b x — a x )/N x , and Ay = (b y — a y )/N y . Therefore, we take as our 
relative L2-norm indicator the following easily computable quantity: 



i=l Z^i = l Z^fe=l 



E 2 (Ax,Ay,M) := 



■>M(M+l)/2 



7 (k) 
ij 



(fc) 
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